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HEATING THE INTRA-CLUSTER MEDIUM BY JET-INFLATED 

BUBBLES 


Shlomi and Noam Soker^ 

ABSTRACT 

We examine the heating of the intra-cluster medium (ICM) of cooling flow clusters of galaxies 
by jet-inflated bubbles and conclude that mixing of hot bubble gas with the ICM is more impor¬ 
tant than turbulent heating and shock heating. We use the PLUTO hydrodynamical code in full 
3D to properly account for the inflation of the bubbles and to the multiple vortices induced by the 
jets and bubbles. The vortices mix some hot shocked jet gas with the ICM. For the parameters 
used by us the mixing process accounts for about four times as much heating as that by the 
kinetic energy in the ICM, namely, turbulence and sound waves. We conclude that turbulent 
heating plays a smaller role than mixing. Heating by shocks is even less efficient. 


INTRODUCTION 


A negative feedback mechanism determines the thermal evolutio n of the intra-cluster medium (IC M) in 


the inner regions of cooling flow (CF) clusters and gro ups of galaxies (|McNamara fc Nnlsenll2QQ7ll2QI2 


some 


examples include Farage et al. I boii Ipfromm^boish . This feedback me chanism is driven by active galactic 


nucleus (AGN) jets that in fl ate X-ray deficient ca vities (bubbles; e.g., Dong et al. 2 QIq[ O’Sullivan et al 


2011 


Caspar! et al 


2052 (iBlanton et ah, __ 

doitti et, alJl^riT . among many others. 


20I2allbl : Gilkis fc Soker 2012 ). _ Exam ples of bubbl es in cooling flows include Abell 

2nnr NGC f)338 (IPa,ndgeet a,l.l 120121) ■ NGC 5044 dPa.vid et al.lboOfll) . and HCG 62 


evolution, their interaction with the ICM, and the dynamics and formation of cold regions (e.g. 

Omma et al. 

2004 

: Heinz & Churazov 

2005: 

Roediger et al. 

2007 

: Sternberg & Soker 

2008b: Gilkis & Soker 

2012|). Some 


cool regions further cool and flow inward to feed the AGN. The process of feeding the AGN with cold 
clumps in the feedback mechanism cycle is termed the cold feedback mechanism^ and was suggested by 


gas and bv more detailed studies (e.g.. Revaz et al. 2008: Pope 2009: Pizzolato & SokerlI 2 OIO: Edge et al. 


2010 

: Wilman et al.lboil: Nesvadba et al.ll20 

III Gavagnolo et al.l[20II:lGaspari et al. 

20I2allb:lMcGourt et al. 

2012 

: Sharma et al.l 

I 2 OI 2 I: lEarage et al. 2012 

: Wagh et al.l 20I4I: Baneriee & Sharma 

boil: McNamara et al. 

2014 

: Voit & Donahuell20I5: Voit et al. 20I5l 

Li et al.ll20I5:IPrasad et al. 20I5l Russell et al. 20I5l Tremblav et al. 

2015 

: lEogartv et al.l 2015). 


Gilkis fc Sokerl (|20I2h and iHillel fc Sokerl (|20I4l ) have shown that the mixing of the hot bubble gas with 
the ICM is a much more efficient heating process than heating by repeated shocks; hereafter termed shocks- 
heating. As the axis of the bubbles changes with time and/or the central j ets’ source and the ICM move 
relative to each other, e.g., the galaxy group NGC 5813 (jRandall et al.ll2QI5[) . over time the mixing-heating 
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operates in all directions. Sh oc ks-heating; is extreme ly inefficient, and is unlikely to explain the feedback 


mechanism (|Soker et alJl2Q13[) . iRandall et alJ (|2 01 51) argue for shocks- heating in NGC 5813, ignoring the 


counter arguments for shocks-heating in NGC 5813 (ISoker et al.N2Q13r). W e hold that even in NGC 5813 
mixing is more important than shocks in heating the ICM ( Soker et al. [2015 1. This heating by mixing closes 


the cycle in the cold feedback mechanism. 

A key to the study of the mixing heating is to inflate bubbles se lf-consistently. This requires either 
slow ( s ub-relativistic) massive wide ( SMW) jets (|Sternberg et al. 1120071 ). precessing let s (ISternberg fc Soker 


2008a 


Falceta-Goncalves et al 


_±J.S _ 

2 OI 0 I 1 . or a relative mot ion of the jets to the medium ( Brhggen et al. 2003 


Sokerl 20^ Morsonv et al. 2010[ Mendvgra! 


et al.N2012l). In the present study we inflate bubbles by SMW 


Moe et al.l 120091: iDnnn et al.l I 2 OIOI: lArav et al.l 120131) . but our 


jets, also supported by observations (e.g., 
results hold for bubbles inflated by precessing jets or a relative motion of the ICM as well. 


In a recent paper I Zhuravleva et al.l (|20I4f ) claim that heating by turbulence dominates in the Perseus 
cooling flow cluster. We here compare the energy that jet-inflated bubbles ch annel to turbul e nce a nd to 
the direct the r mal e nergy of the ICM. The axi-symmetric 2D simulations of ICilkis fc Sokerl (|20I2[) and 
Hillel fc Sokerl (|20I4l ) are insufficient to account for the vortices and turbulent motion, which are expected 
to be fully three-dimensional. Therefore, as 2D axi-symmetry restricts the motion and changes the character 
of the turbulence invoked, for the goal of this paper, that is to study the role of turbulent heating, we lift 
this restriction and employ 3D simulations. 

The 3D numerical setting is described in section [2l In sections [3] we present the flow structure, and in 
section m we study the heating of the ICM. Our short summary is in section [5l 


2. NUMERICAL SETUP 


We use the PLUTO code ( Mignone et al. 200?!) for the hydrodynamic simulations, in a three-dimensional 
Cartesian grid with adaptive mesh refinement (AMR). The computational grid is in the octant where the 
three coordinates x, y and 2 ; are positive. At the x = 0^ y = 0 and z = 0 planes we apply reflective boundary 
conditions. The ^ coordinate is chosen along the initial axis of the jets. In reality two opposite jets are 
launched simultaneously, such that the flow here is assumed to be symmetric with respect to the 2 : = 0 
plane, amounting to reflective boundary conditions at 2 ; = 0. The base computational grid (lowest AMR 
level) spans the cube 0 < x,^, 2 ; < 50 kpc, with 16 divisions in each direction. Up to 5 AMR levels are 
employed with a refinement ratio of 2. Thus, the highest resolution is ~ O.I kpc. The refinement criterion 


is the default AMR criterion in PLUTO v. 4.1, based on the second derivative error norm (|Lohnerl 119871) of 
the total energy density. 

Heat conduction and viscosity are not in cluded in th e simulations. However, local heat conduction is 
expected to be efficient on scales of < O.I kpc ( Soker 2 QIq[) . and therefore is effectively incorporated via the 
resolution of the simulations and does not need to be included. The conclusions are found not to be sensitive 
to the resolution. Magnetic fields, also not included in the current numerical study, would in reality prevent 
global heat conduction, but not local heat conduction. The vortices would entangle field lines, leading to 
local magnetic field line reconnection and hence allowing local heat conduction. 

In this study we examine the increase in the kinetic energy of the ICM. The kinetic energy results from 
turbulence, sound waves, and gas mo tion induced by cavity m otion, i.e., cavity heatin g dNulsen et al. 2007 ). 
The dissipation of turbulence (e.g., IZhnravleva et al.N2QI4l) and sound waves (e.g., IShabala fc Alexander 
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20071) heats the gas. We do not include dissipation, but we calculate the increase in the kinetic energy of 
the ICM and term it in short turbulence. However, it should be kept in mind that it actually includes the 
energy of sound waves and uplifted gas, hence we do account for sound wave heating and cavity heating. If 
the dissipation is not 100% efficient in transferring kinetic to thermal energy, then we are overestimating this 
turbulent energy source for heating. The heating by shock wave s is examin e d dire ctly, as shock dissipation 
is treated by the numerical code. Heating by cosmic rays, e.g., iGno fc OhI (|2008l ). are not included in our 
simulations at all. 

On the outer boundaries we use an outflow boun dary condition. At the boundary ^ = 0 we inject into 


the grid a jet with a half-opening angle of 0-^ = 70° (|Sternberg et al. 1120071) . The jet material is inserted 


through a circle x‘^ + y‘^ < r? on the z = 0 plane with a radius of rj = 3 kpc. The initial jet velocity in 
the nominal case is v-^ = 8200 km s“^, a Mach number of about 10. The direction of the velocity at each 
injection point (x, y, 0) in the circle is h = (x, y^h^)/y‘^ -\- h?, where hj = r-J tan^j. The jet is injected 
during each active episode, and when the jet is turned off reflective boundary conditions apply in the whole 
z = 0 plane. In the nominal run, the jet is periodically turned on for 10 Myr and off for 10 Myr. The power 
of the two jets together is 

E 2 i=‘ 2 x 10^® erg s'S (1) 


half of it in each direction. The mass deposition rate is thus 

2 P' 

M2j = ^=94Mo yr-i. 


( 2 ) 


In two additional simulations we vary two parameters. In one simulation, termed Run B below, we 
reduce the jet activity and quiescence time intervals by a factor of 5. In a different simulation. Run C, we 
increase the velocity by a factor of and lower the mass deposition rate, keeping the jet power unchanged. 

As we do not include magnetic fields, hence their influence on instabilities, the roles of magnetic fields in 
the mixing process and in small scale motion s are not explored he re. Although in two dimensions magnetic 
fields act to stabilize some instabilities fe.g.. iRobinson et al.ll2QQ4l) . in three dimensions the magnetic fields 
cannot suppress instabilities with wave-vector components perpendicular to the field lines. Hence, we do not 
expect magnetic fields to suppress mixing, but only to influence the details of the mixing. Another effect is 
that the violent vorticity we find can amplify magnetic fields to the point that they start to reconnect and 
heat the gas. This overall process channels kinetic energy to magnetic energy and then to thermal energy, 
hence increasing the heating efficiency of the mixing process. Overall we do not expect magnetic fields to 
influence the main conclusions of this study. 

The simulation begins wi th an isothermal box of gas at an initial temperature of Ticm(O) = 3 x 10^ K 
with a density profile of (e.g., Vernaleo fc RevnoldsI 2006 ) 


picuir) = - 


Po 


1 + {r/ay 


-|3/4’ 


(3) 


with a = 100 kpc and po = 10“^^ g cm“'^. A gravity field is added to maintain an initial hydrostatic 
equilibrium, and is kept constant in time. We include radiative cooling i n the simulations, where the t abulated 
cooling function is taken from the solar-metalicity values of Table 6 in [Sutherland fc Dopital (| 19931) . 


We run simulations with either 4 or 5 levels of AMR, to check for numerical convergence. We found 
no noticeable differences in the results. We also run simulations with a twice as large grid and no reflecting 
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boundary conditions on the plane x = 0. Namely, the x-coordinate covers the range —50 kpc < x < 50 kpc. 
We found no significant changes in the results, hence conclude that the reflecting boundary conditions do 
not affect our results, both qualitatively and quantitatively. 


3. FLOW STRUCTURE 

We begin by presenting in Fig. [T]the flow structure at different times. This high-resolution simulation 
was run for 50 Myr. To follow the evolution of the same setting for 240 Myr we used a lower resolution grid. 
In the common time the results of the two runs are similar. 
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Fig. 1.— Evolution of the flow for the nominal case parameters given in section[2l presented in the meridional 
plane ^ = 0 at four times. The jet is injected through a circle of radius rj = 3 kpc on the 2 ; = 0 plane. The 
half-opening angle is 6>j = 70°, i.e., the jet velocity on the boundary of the circle has an angle of 70° with 
respect to the 2 ;-axis. The left panels present temperature maps and the right panels mass density maps. 
The color coding of temperature and density is in logarithmic scale. Arrows show the velocity, with length 
proportional to the velocity magnitude. A length of 1 kpc on the map corresponds to 1700 km s“^. When 
the jet is active, the length of arrows close to the origin corresponds to 8200 km s“^. 
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As evident from Fig. [TJ e.g., at t = 32 Myr, we manage to inflate bubbles which match those seen i n 


observations as X-ray deficient cavities (|Sternberg et al. 1120071: iHillel fc Sokerll2Ql4 iGilkis fc Sokerll2Q12[) . 


The inflation of bubbles similar to those observed in a self-consistent manner is crucial to the study of the 
different heating mechanisms: mixing, shocks, and the excitation of ICM turbulence that can dissipate later 


In order to compare to observational data, in the left panel of Fig.[2]we present mock X-ray images of the 
nominal simulation at time t = 50 Myr. This image was created by folding the simulated octant twice, such 
that the four quarters of the plane are identical. Although the initial setup of the simulation is cylindrically 
symmetric, hydrodynamic instabilities magnified numerically cause departures from that symmetry. In 
particular, the surface of the low density hot bubble is ripply and uneven, with several relatively dense 
filaments crossing that volume near the bubble’s surface (later discussed in relation to Fig. [ 7 ]). These appear 
as a net of bright filaments covering the dark (low density) bubble on the left panel of Fig. [2j Present X-ray 
observations are not sharp enough to detect such filaments. In addition, more sophisticated simulations 
might have erased these filaments by including realistic local heat conduction. 


I 
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Fig. 2.— Mock X-ray image (left panel) and emission-weighted temperature map (right panel) of the nominal 
simulation at time t = 50 Myr. The simulated octant was folded twice to create these panels. Hence the 
four quarters about the origin are identical to each other. The X-ray image is created by the integration 
of the density squared along lines parallel to the ^-axis. The units are arbitrary. The emission-weighted 
temperature map is made by calculating the average, along the line of sight, of the temperature weighted by 
the emission power, Eq. Note the X-ray deficient bubbles inflated by the jet. 


The ‘waist’ between the two bubbles in the plane z = 0 obtained here is wider than in typical observed 
bubbles. Our method of injecting the jet, through a circle on the surface ^ = 0, is one of the causes of the 
wide waist near z = 0. In addition, on the plane 2 : = 0 we have reflective boundary conditions which add 
to the numerical error near the plane. However, overall, we obtain an image which reasonably resembles 
observations. 

The right panel of Fig. [2] is an emission-weight temperature map at the same time, t = 50 Myr. The 
emission-weighted temperature at each point (x, z) in the image is the average of the temperature, along the 
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line of sight (^-axis here), weighted by the emission power, 

_ /A(T)n2Td2/ 
/A(T)n2d2/ • 


(4) 


For reference we also provide representative shock Mach number values. At the front of the jet, on the 
z-axis, the first shock wave running through the ICM has a Mach number of = 1.53 at z = 20 kpc and 
a Mach number of A4 = 1.18 when it reaches a distance of 2 : = 40 kpc. As we show later, such weak shocks 
are insignificant in heating the ICM. 


To better capture some properties of the flow we present artificial flow quantities called ‘tracers’. The 
tracers are frozen-in to the flow, and hence represent the spreading with time of gas starting in a certain 
volume. A tracer’s initial value is set to ^(0) = 1 in a certain volume and ^(0) =0 elsewhere. If the traced 
gas is mixed with the ICM or the jet’s material, its value drops to 0 < ^(t) < 1. 


In Fig. [ 3 ] we present maps of a tracer that follows the ejected jet material. In the simulation the jet is 
injected via the boundary condition at ^ = 0 (see sectionj^), and the injected material is marked with a tracer 
that follows its evolution throughout the simulation. We plot the tracer distribution in the ^ = 0 meridional 
plane at two times as indicated in the different panels. Each point in space in which the tracer value is 
^ > 0 indicates that material from the jet has reached it. We And again that a prominent phenomenon in 
such simulations of self-co nsistent bubble inflation is the appearance of vortices and mixing (|Gilkis fc Soker 


2012 


Hillel fc Sokerll2QI4l) . The inflation of the bub ble by the jet induces a complicated flow stru cture with 


multiple vortices, some induced by vortex shedding (|Refaelovich fc Sokerll2QI2l: IWalg et al.ll2QI3[) . 


In the upper-left left panel of Fig. [3] we present the entire range of 0 < (^ < I, while in the upper-right 
panel we zoom on a specific region and emphasize the range 0 < ^ < O.I. The arrows represent the direction 
and magnitude of velocity. Too long arrows are truncated with a length corresponding to a maximum velocity 
Vm as explained the caption. Vigorous mixing caused by vortices can be seen in these two panels, i.e., the 
vortex at (x, z) = (5,18) kpc. In the upper-right panel we can notice jet’s material that reaches the point 
(x, z) = (23, 3.5) kpc. This clearly shows that the mixing of hot shocked jet’s gas can heat the ICM to large 
distances in directions perpendicular to the initial jet direction. Considering that in a more realistic scenario 
the jets precess and eject material in other directions, the mixing is more efficient even than what is found 
here. 


The two lower panels of Fig. [3] present a later evolution time. The curved arrows in the lower-left panel 
represent stream lines, and the arrows in the lower-right panel represent mass flux (j) = vp. Both the stream 
lines and the arrows show a shock expanding outward at r ~ 35 kpc. We can see again jet’s material that 
reaches regions away from the jets’ axis, here at (x, z) = (28, 20) kpc. 

In Figs, m and [ 5 ] we show the evolution of tracers that are frozen-in to the ICM. In both cases the tracer 
follows the gas that started in a torus around the ^ axis, and the radius of the cross section of the torus is 
2.5 kpc. In Figs. |4]we follow the tracers whose torus cross section is centered at (x, z) = (10, 5) kpc, and in 
Fig. [5] whose torus cross section is centered at (x, z) = (20,15) kpc. The bottom right panel of Fig. [5] shows 
the tracer of the jet rather than the tracer of the ICM. Both figures clearly show the efficient mixing of the 
ICM with the hot bubble gas. 

The tracer in Fig. 0] is initially pushed outward by the outgoing shock waves and the growing bubble of 
hot gas that induces sound waves. However, at t ^ 40 Myr it begins to be mixed with hot jet gas, as the 
large vortex induced by bubble inflation tears it apart. At t ^ 100 — 150 Myr it becomes almost completely 
mixed with its surroundings as well as with the jet gas. 
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Fig. 3.— A ^ = 0 slice of the tracer that is frozen-in to the jet material, at times t = 50 and 100 Myr in 
a simulation with the same parameters as in Fig. [H The value of the tracer is color-coded according to the 
scale shown on the right of each panel. In two panels the length of the arrow is linear with velocity up to 
an upper limit In the upper-left panel the largest velocity arrow corresponds to = 8200 km s“^, the 
initial jet velocity. Faster regions are presented with an arrow length corresponding to = 8200 km s“^. 
The upper-right panel zooms on a region away from the center to emphasize mixing of jet’s material away 
from the center and near the equatorial plane. The largest velocity arrow corresponds to = 3000 km s“^. 
In the lower-left panel we show stream lines of the flow at t = 100 Myr, and in the lower-right panel we show 
mass flux arrows, i.e., (j) = vp. A length of 1 kpc on the map corresponds to 0 = 2.3 x 10“^^ km s“^ g cm“^. 

The tracer of the gas starting further out and presented in Fig. 0 follows a similar pattern, but begins to 
be mixed with the hot bubble gas at much later times. Even that it is pushed out to ~ 30 kpc, it eventually 
flows inward and mixes with the hot bubble gas and heats up, as we show in section 01 In the lower-right 
panel of Fig. [5] we show the tracer of the jet. The two lower panels that are given at the same time clearly 
present the efficient mixing of the ICM and the hot bubble gas. For example, a clear vortex mixing the two 
media is seen circulating around (x, z) = (17, 5) kpc. 

An inflow of gas near the equatorial plane is seen at the latest time in Fig. [5l This carries ICM gas 
















































































































- 8 - 




15 20 25 

X (kpc) 



1.0 

0.9 

0.8 

0.7 

0.6j 

CO 

0.5 2 

CD 

O 

0 . 4 ^ 

0.3 

0.2 

0.1 

0.0 



0 5 10 15 20 25 30 35 40 

X (kpc) 


Fig. 4.— Evolution with time of the gas that at t = 0 was contained in a torus whose axis is the 2 ;-axis 
and whose cross section is a circle centered at (x, z) = (10, 5) kpc with a radius of r = 2.5 kpc. Shown 
is the concentration of this gas, as followed by its tracer, in the meridional plane ^ = 0 at four times. 
Flow parameters are as in previous figures. The largest velocity vector corresponds to = 400 km s~^ 
0.5A1icm, where Micm is the Mach number in the ICM. Higher velocities are marked with arrows with the 
same length as that of Um- 


toward the jets, followed by mixing and heating. This demonstrates again that the jets can heat material 
even in regions further out and near the equatorial plane. If the jets cont inue to be active fo r a very long time 
along the same axis, a massive meridional flow will develop in the ICM (|Soker et al.ll2015[) . This flow carries 
cooler ICM gas towards the center and mixes it with the hot bubble gas. This is an indirect mechanism by 
which a jet may heat up gas in regions further away from its axis. 


4. HEATING THE ICM 

We turn to present the energy history of the ICM. In order to follow specific regions we mark them 
with ‘tracers’ as explained in section [3l The total mass of a traced gas is given by the sum Mtr = 
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over all grid cells. As the tracers are advected with their associated mass, the total traced mass is constant 
with time, as long as traced gas does not leave the grid. Since there are outflow boundary conditions on the 
grid faces x = 50 kpc, ^ = 50 kpc and z = 50 kpc, and traced gas is dragged by the jets after some time, 
traced gas does eventually leave the grid. When that happens we stop following the thermal evolution of 
that traced gas. 

We use the tracers to define the total energy Et^ of a traced region as 

Etr = E,i^iEi, (5) 

where i runs over all grid cells. In this work, Etr stands for the traced gas’s total kinetic energy, internal 
energy, or gravitational energy, and Ei stands for the corresponding energy in cell i. In Fig. [6] we present 
the energy histories of different traced regions in the simulations. The internal energy and the gravitational 
energy are drawn with respect to their initial value which is different than zero. The initial kinetic energy is 
zero. The initial thermal energy of each traced region is given in the caption. 

In addition to the nominal jet parameters specified in sections [2] and [3l we conducted two simulations 
with varying parameters, in order to check the robustness of our conclusions. The nominal simulation is 
represented in the three panels in the left-hand column of Fig. [6l Two additional simulations. Run B and 
Run C, are represented in the middle and right-hand columns respectively. We shall elaborate on them 
below. 

Consider first the results of the nominal simulation, i.e., the left-hand column. In the panels in the top 
row of Fig. [6] we present the energy evolution of mass located at t = 0 inside a sphere of radius r = 15 kpc 
centered at the origin. Because of the structure of our grid, only an eighth of the sphere is traced. The 
initial internal energy of this traced gas in the simulated eighth is F^in(O) = 3.1 x 10^^ erg, such that the 
maximum value of the extra thermal energy in the plot is A = [F^in(max) — Fin (0)]/Fin(0) = 390% (in the 
nominal case) above the initial internal energy. The calculation is stopped when traced gas starts leaving 
the computational grid. 

The panels in the middle row present the energy history of the traced gas whose spreading history is 
displayed in Fig. 01 In this case Fin(O) = 5.5 x 10^^ erg, and A = [£’in(max) — Fin(0)]/Fin(0) = 550% (in 
the nominal case). The calculation is stopped when traced gas starts leaving the computational grid. The 
panels in the bottom row present the energy history of the traced mass whose spreading history is shown 
in Fig. [5l In this case Fin(O) = 1.1 x 10^^ erg, and A = [Fin (max) — Ein(0)]/Fin(0) = 50% (in the nominal 
case). 

A general characteristic behavior emerges from the energy evolution of the traced masses presented here, 
as well as of other traced regions we have studied but left out of this text. The variation on a time scale of 
20 Myr corresponds to the time period of the jet activity. Fach active phase lasts for 10 Myr, followed by a 
quiescence phase of 10 Myr. The periodic jet activity induces gas motion in the ICM such that the kinetic 
energy of the ICM at each point varies with time. This variation is more or less in a periodic manner, with 
only a moderate increase in the kinetic energy over several periods. The gravitational energy changes by a 
negligible amount compared to both the kinetic and internal energies. The main energy increase over several 
periods comes from an increase of the internal energy. 

Consider the middle panel in the left-hand column of Fig. [6l which shows the energy history of the 
tracer presented in Fig. 0) There are some variations due to two shock waves that pass through the traced 
region. The shocks increase the internal energy, but the gas re-expands and the net thermal energy gain 
is negligible. This traced mass begins to mix with hot bubble gas at t ~ 40 Myr. From that time on the 
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mean internal energy increases significantly. The upper right panel of Fig. 0 ] shows the traced gas at time 
t = 50 Myr, while both upper panels of Fig. [3] show the presence of the jet material at the same time. It is 
apparent that the hot shocked jet (bubble) gas is being mixed vigorously with the ICM gas marked by the 
tracer. The mixing is facilitated by the vortices induced by the jet-infiated bubble. 

Similar heating by mixing happens in the tracer presented in Fig.O whose energy history is shown in the 
bottom panel in the left-hand column of Fig. [6l Periodic heating by shocks is clearly seen, although over a 
long time shocks cannot compete with radiative cooling. Indeed, at early times the internal energy decreases. 
Here the mixing begins at t ^ 110 Myr. In the t = 100 Myr panel of Fig. Owe can see a vortex that is just 
about to touch the traced region and mix it with hot bubble gas that can be seen in the t = 100 Myr panel 
of Fig. O 

The heating is mainly due to mixing of the ICM gas with shocked jet material, i.e., the hot bubble 
gas. Both the mixing that leads to an increase in internal energy and the turbulent energy are driven 
by vortices. It is hard to estimate directly the amount of turbulent kinetic energy of the ICM from the 
simulations, because the grid resolution is not high enough to capture all relevant scales. However, we give 
rough arguments that the kinetic energy we calculate is not much smaller than the corresponding True’ 
turbulent energy. Firstly, the kinetic energy is calculated with respect to the center of mass of the whole 
system, and not with respect to a local mean flow of the gas. Thus, it overestimates the turbulent energy 
as it is usually defined. Secondly, our calculated kinetic energy does not properly include the kinetic energy 
in smaller scales than the grid resolution, and thus it underestimates the turbulent kinetic energy. However, 
the kinetic energy is mainly due to vortices which are significantly larger than the grid resolution. Therefore, 
we do not expect that the underestimation of the turbulent kinetic energy in small unresolved scales is 
significant. Thirdly, some of the kinetic energy we count is due to sound waves that propagate out and do 
not heat the inner region. Finally, motion in smaller scales is less energetic, and it also tends to dissipate to 
internal energy. 

We now turn to the additional two simulations presented in Fig. [6l In one simulation. Run B, shown 
in the panels in the middle column of Fig. [6l the jet is periodically turned on for 2 Myr and off for 2 Myr, 
instead of 10 Myr/10 Myr in the nominal run. Density, temperature, and velocity maps in the ^ = 0 plane 
at time t = 50 Myr of Run B are shown in the middle column of Fig. 0 In Run C, shown in the panels in 
the right-hand column of Fig. [6l the jet periodicity is unchanged, but the mass deposition rate is decreased 
by a factor of 10 and the jet velocity is increased by a factor of with respect to the nominal simulation. 
The rest of the parameters are as in the nominal simulation. The flow maps for Run C in the ^ = 0 plane 
at t = 50 Myr are shown in the right-hand side column of Fig. 0 In all simulations the mean jet power is 
the same, Eq. O- 

The detailed flow and vortex structure in the three simulations is dependent upon the parameters in 
a complicated way. However, we can draw conclusions which do not depend on the exact value of the 
parameters. In Run B (middle column of Fig. [6]), the time period of jet episode activity is 5 times shorter. 
The immediate effect seen is the smaller and more frequent ‘bumps’ as a result of the shock waves. In 
addition, the shape of the hot bubble is elongated in the direction of the axis of the jet, while maintaining a 
similar volume which depends only on the jet energy output. Run C, with the higher jet velocity and lower 
mass deposition rate, has a similar general structure to the nominal run, but the details are different. In Run 
C there is less momentum per unit mass and per unit energy, compared with the nominal run. The bubble 
front along the 2 ;-axis after 50 Myr extends only to 27 kpc, compared with 32 kpc in the nominal run. The 
vortex on the right side of the bubble in Fig. 0 is larger in Run C than in the nominal run. The waist in Run 
C is wider than in the nominal run. The nominal run better matches the morphology of observed bubbles. 
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This is one of the reasons we u se slow lets (ISternberg et al. 1120071) . Another reason is recent observations of 


such outflows from AGN (e.g., Chamberlain et al. 20151) 


The energy histories of the traced regions in these simulations, Fig.[6j is similar to the nominal run with 
the following exceptions. In Run B (middle column of Fig. [6]), in the tracer shown in the bottom panel, the 
traced gas gets close to the hot bubble at t = 100 — 150 Myr but does not mix with the hot gas. The more 
frequent episodes create less vortices and in fact push out the traced gas almost continuously. In contrast, 
in the nominal run, where the jet is periodically turned on and off for 10 Myr, the longer quiescence time 
allows for the cooled ICM gas to be pulled in to the center and thus to heat up by mixing with the hot 
bubble gas. In Run B this happens to a lesser extent. Thus, the traced ICM gas remains in the dense shell 
surrounding the hot bubble, which in turn radiates its energy more efficiently. This radiative cooling is the 
cause of the decline of internal energy starting at time t ~ 180 Myr. As seen in this example, more shock 
waves (keeping the jet power fixed) may, in some cases, actually heat some of the ICM gas less efficiently, 
because it pushes the gas out more continuously, which inhibits mixing. 


Similar cooling of the same traced region happens in Run C, as seen in the bottom-right panel of Fig. [6l 
Here, the gas begins to mix at time t ~ 110 Myr, as in the nominal simulation. However, the slightly 
different vortex structure does not pull in the entirety of the traced gas and most of it remains out of the 
hot bubble, in the surrounding dense shell. Thus, as in Run B, it loses internal energy quickly by radiative 
cooling. This may be a mechanism for feeding the AGN by cold gas. ICM gas gets pulled in to the center 
during quiescence periods because of the lower pressure in the hot bubble, cools radiatively when its density 
increases, and is then pulled in further. At this stage it either heats up or cools further and falls in to feed 
the AGN. 


The main results from the comparison of the different cases is as follows. Because of the different flow 
structure in simulations with varying parameters, some of the gas heats up more and some less, depending 
on the location of the gas, but mixing is the main heating process among those probed in the simulations. 

We thus conclude by arguing that mixing, due to vortices, is more important as a heating process of the 
ICM gas by jets than turbulent heating and shock heating. Shock waves induce an expansion and contraction 
and increase kinetic energy by creating motion, but their contribution to the heating is largely temporary 
and almost periodic in our setting. Turbulent motion is also not the main heating mechanism, although it 
might carry up to ~ 0.2 of the energy transferred from the jets to the ICM. 


5. SUMMARY 


Motivated by the recent claim of I Zhuravleva et al.l (|2Q14l) that turbulent-heating can counter radiative 
cooling in the cooling flow clusters Virgo and Perseus, we compare the energy channeled from jets to ICM 
turbulence with the direct thermal heating by mixing hot bubble gas with the ICM. We have done so by con¬ 


ducti ng 3D hydrodynamical simulations of wide jets that inflate bubbles. We used the PLUTO (|Mignone et al 


20071 ) hydrodynamical code in Cartesian coordinates. 


The inflation of bubbles leads to the development of vortices in the entire volume around the expanding 
jets and bubbles, as can be see in Figs. [T]l5] and [71 These vortices excite turbulence in the ICM. These 
vortices also efficiently mix ICM gas with the hot bubble gas. This mixing can best be seen in Fig. 0] where 
the tracer (colored red) of ICM gas located initially within a torus is seen to be mixed with the bubble gas. 

We use slow massive wide (SMW) outflows. Such flows are supported by observational findings (e.g.. 
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Chamberlain et alJ l2Q15l and references therein). In many cases outflows from AGN include relativistic 
components. By comparing the results of Run C with the nominal run, e.g., in Fig. [71 we note that when 
the velocity is three times as high, the bubble stays wide and a large vortex exists on its side. Higher 
velocities do not change the mixing part by much when wide (fat) bubbles are inflated. Therefore, we do not 
expect that relativist outflow will change our present conclusions, as long as wide bubbles are inflated. For 
that, either the jets are wide, or there is a relative transverse v elocity between the jets a nd the ICM, e.g., 
motion of the AGN relative to the ICM and/or jets’ precession (|Sternberg fc Sokerll2QQ8al) . Well-collimated 
relativistic jets without transverse velocity relative to the ICM will not heat the ICM efficiently away from 
their propagation cone. They will penetrate to large distances. They can form radio lobes at large distances, 
but will not heat the ICM close to the center. 

We then compared the increase in thermal energy of the ICM medium we trace, with the increase of its 
kinetic energy. We found that the energy channeled to directly heating the gas is more than 3 times larger 
than that channeled to kinetic energy of the ICM, and typically larger than 4 times the kinetic energy. This 
is our main result and it is presented in Fig. O As only part of the kinetic energy will turn to turbulence 
(some is just large-scale motion), we conclude that turbulent-heat ing is < 30% as efficient as mixing-heating . 
Heating by shocks is very small, as was shown in 2D simulations (|Gilkis fc Sokerll2QI2l:lHillel fc Sokerll2QI4[) . 
and was reinforced here in 3D simulations. The significance of the new 3D simulations is that the vortices and 
turbulent motion is not artificially confined to be cylindrically symmetric. Although the setup is cylindrically 
symmetric, small perturbations and numerical errors are magnified by hydrodynamic instabilities, and the 
flow is fully three-dimensional. This is the main difference observed in the 2D and 3D simulations. The 
conclusions regarding the heating mechanisms, nonetheless, are similar. 

It is important to note that since by our finding turbulence carries < 25% of the energy deposited 
by the jets in the inner region, turbulence may s till be non-neg l igible in the IGM, as was found recently by 
Zhuravleva et al. (2014), Zhuravleva et al. (2015), Arevalo et al. ( 2015h . and marginally by Anderson fc Sunvaev 
(|2015h . However, since only a portion of the kinetic energy develops to turbulence, and since the total injected 
jet energy can be channeled into other forms of energy, in some cases turbulence might be insignificant. 


Zhuravleva et al.l (|2014l) find observationally that turbulence is present and significant in the ICM, and 
that the rate of turbulent heating matches the rate of radiative cooling. In contrast, in our simulations 
we have found that heating by mixing is more significant. Our results further show that a non-negligible 
amount of energy ends up as kinetic energy of the ICM. However, substantial heating occurs only when 
mixing takes place. In addition, we find large density fluctuations induced in the ICM by the mixing process 
and bubble motion; this is best seen in the upper panels of Fig. [71 It seems that a large fraction of the 
kinetic energy is not turbulent motion that dissipates. The kinetic energy and density fluctuations are 
associated also with the bubble acti vity and global gas motio n; only p art of of the kinetic ener gy belongs to 
the universal turbulent cascade that [Zhuravleva et al.l (|2QI4l ) refer to. [Zhuravleva et al.l (|2QI4f) write: “It is 
difficult to prove unambiguously that we are dealing with a universal turbulent cascade, as other structures 
(e.g., edges of the bubbles, entrainment of hot bubble matter, sound waves, mergers and gas sloshing) might 
also con tribute to the observed density-fluctuation spectra.” We suggest that the discrepency between the 
claim of [Zhuravleva et al.l (|2QI4l) for efficient turbulent heating and our finding results from a large fraction 
of the density perturbations and turbulent energy that is indeed in the other structures listed by them. This 
definitely deserves further investigation. 


We thank an anonymous referee for very valuable suggestions and comments. 
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Fig. 5.— Evolution with time of the gas that at t = 0 was contained in a torus whose axis is the z-axis and 
whose cross section is a circle centered at (x, z) = (20,15) kpc with a radius of r = 2.5 kpc. In the first five 
panels the concentration of this gas, as followed by its tracer, is shown in the meridional plane ^ = 0 at five 
times. The lower-right panel shows the tracer of the jet. Flow parameters as in previous figures. The largest 
velocity vector corresponds to Um = 400 km s“^, a Mach number of about 0.5. Higher velocities are marked 
with arrows with the same length as that of Um- 
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Fig. 6.— The energy history of different traced regions of the nominal simulation studied here as well as in 
two additional simulations. The three panels in the left-hand column represent the nominal simulation. In 
the second simulation, Run B, presented in the middle column, the jet is periodically turned on for 2 Myr 
and off for 2 Myr. The rest of the parameters are as in the nominal simulation used in the rest of this work. 
In the third simulation. Run C, presented in the right-hand column, the jet mass deposition rate is decreased 
by a factor of 10 and the jet velocity is increased by a factor of ^/TO with respect to the nominal simulation, 
such that the jet’s power is unchanged. Otherwise, the setup is the same. In the three panels in the top row 
we present the evolution with time of three energy components of the ICM gas that starts (before the jets 
become active) inside an eighth of ball with r = 15 kpc centered at the origin. It is an eighth of a ball as 
we simulate one eighth of the space. The green line in each panel represents the internal energy, the blue 
line represents the kinetic energy, and the red line represents the gravitational energy of this traced gas. 
The middle-row panels show the energy histories of the torus shown in Fig. 01 and the panels in the bottom 
row show the energy histories of the torus shown in Fig. [5l All energies are shown relative to their values 
at t = 0. The initial internal energies F^in(O) of the traced regions, from top to bottom, are 3.1 x 10^^ erg, 
5.5 X 10^^ erg and 1.1 x 10^^ erg, respectively. The top- and middle-row panels are cut off at the time when 
the traced material starts leaving the grid. 
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Fig. 7.— Mass density and temperature maps of the nominal run (left column), Run B (middle column) 
and Run C (right column) whose energy histories are presented in Fig. O The ^ = 0 slices of the simulations 
are presented at time t = 50 Myr. The color coding of the density is linear and presented in the range 
(0.5,1.5) X 10 g/ cm^, and the temperature coding is in logarithmic scale. Arrows represent the velocity, 
where a length of 1 kpc on the map corresponds to 1700 km s“^. 










